Take_home_Ex1

Modified

November 26, 2023

Overview

Objectives

Load Packages and Data

Load packages

pacman::p_load(sf,tmap,spdep,sf,tidyverse, dplyr, mapview)

Load data

Import aspiatial data

odbus <-read_csv("data/aspatial/origin_destination_bus_202310.csv")

Getting the data of the different duration

weekdayAM_6_9<-odbus %>% 
  filter(DAY_TYPE == "WEEKDAY") %>%
  filter(TIME_PER_HOUR >= 6 &
           TIME_PER_HOUR <= 9) %>%
  group_by(ORIGIN_PT_CODE) %>%
  summarise(TRIPS = sum(TOTAL_TRIPS))
weekdayPM_5_8<-odbus %>% 
  filter(DAY_TYPE == "WEEKDAY") %>%
  filter(TIME_PER_HOUR >= 17 &
           TIME_PER_HOUR <= 20) %>%
  group_by(ORIGIN_PT_CODE) %>%
  summarise(TRIPS = sum(TOTAL_TRIPS))
weekendAM_11_14<-odbus %>% 
  filter(DAY_TYPE == "WEEKENDS/HOLIDAY") %>%
  filter(TIME_PER_HOUR >= 11 &
           TIME_PER_HOUR <= 14) %>%
  group_by(ORIGIN_PT_CODE) %>%
  summarise(TRIPS = sum(TOTAL_TRIPS))
weekendPM_4_7<-odbus %>% 
  filter(DAY_TYPE == "WEEKENDS/HOLIDAY") %>%
  filter(TIME_PER_HOUR >= 16 &
           TIME_PER_HOUR <= 19) %>%
  group_by(ORIGIN_PT_CODE) %>%
  summarise(TRIPS = sum(TOTAL_TRIPS))

Import geospatial

busstop <- st_read(dsn = "data/geospatial",
                   layer = "BusStop") %>%
  st_transform(crs = 3414)
Reading layer `BusStop' from data source 
  `C:\Users\Lian Khye\Desktop\MITB\Geospatial\geartooth\ISSS624\Take_home_Ex1\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 5161 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 3970.122 ymin: 26482.1 xmax: 48284.56 ymax: 52983.82
Projected CRS: SVY21
glimpse(busstop)
Rows: 5,161
Columns: 4
$ BUS_STOP_N <chr> "22069", "32071", "44331", "96081", "11561", "66191", "2338…
$ BUS_ROOF_N <chr> "B06", "B23", "B01", "B05", "B05", "B03", "B02A", "B02", "B…
$ LOC_DESC   <chr> "OPP CEVA LOGISTICS", "AFT TRACK 13", "BLK 239", "GRACE IND…
$ geometry   <POINT [m]> POINT (13576.31 32883.65), POINT (13228.59 44206.38),…
mpsz <- st_read(dsn = "data/geospatial",
                   layer = "MPSZ-2019") %>%
  st_transform(crs = 3414)
Reading layer `MPSZ-2019' from data source 
  `C:\Users\Lian Khye\Desktop\MITB\Geospatial\geartooth\ISSS624\Take_home_Ex1\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 332 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
Geodetic CRS:  WGS 84
busstop_mpsz <- st_intersection(busstop, mpsz) %>%
  select(BUS_STOP_N, SUBZONE_C) %>%
  st_drop_geometry()
AM_6_9 <- left_join(weekdayAM_6_9 , busstop_mpsz,
            by = c("ORIGIN_PT_CODE" = "BUS_STOP_N")) %>%
  rename(ORIGIN_BS = ORIGIN_PT_CODE,
         ORIGIN_SZ = SUBZONE_C)
PM_5_8 <- left_join(weekdayPM_5_8 , busstop_mpsz,
            by = c("ORIGIN_PT_CODE" = "BUS_STOP_N")) %>%
  rename(ORIGIN_BS = ORIGIN_PT_CODE,
         ORIGIN_SZ = SUBZONE_C)

AM_11_2<- left_join(weekendAM_11_14 , busstop_mpsz,
            by = c("ORIGIN_PT_CODE" = "BUS_STOP_N")) %>%
  rename(ORIGIN_BS = ORIGIN_PT_CODE,
         ORIGIN_SZ = SUBZONE_C)

PM_4_7<- left_join(weekendPM_4_7 , busstop_mpsz,
            by = c("ORIGIN_PT_CODE" = "BUS_STOP_N")) %>%
  rename(ORIGIN_BS = ORIGIN_PT_CODE,
         ORIGIN_SZ = SUBZONE_C)
duplicate1 <- AM_6_9 %>%
  group_by_all() %>%
  filter(n()>1) %>%
  ungroup()

duplicate2 <- PM_5_8 %>%
  group_by_all() %>%
  filter(n()>1) %>%
  ungroup()

duplicate3 <- AM_11_2 %>%
  group_by_all() %>%
  filter(n()>1) %>%
  ungroup()

duplicate4 <- PM_5_8 %>%
  group_by_all() %>%
  filter(n()>1) %>%
  ungroup()
c(duplicate1,duplicate2,duplicate3,duplicate4)
$ORIGIN_BS
 [1] "11009" "11009" "22501" "22501" "43709" "43709" "47201" "47201" "51071"
[10] "51071" "53041" "53041" "62251" "62251" "67421" "67421" "68091" "68091"
[19] "77329" "77329" "82221" "82221" "96319" "96319" "97079" "97079"

$TRIPS
 [1] 17311 17311 11384 11384  1329  1329 25910 25910  6698  6698  4462  4462
[13]  8605  8605   258   258   254   254   337   337  3998  3998   675   675
[25]  3822  3822

$ORIGIN_SZ
 [1] "QTSZ01" "QTSZ01" "JWSZ09" "JWSZ09" "BKSZ07" "BKSZ07" "WDSZ07" "WDSZ07"
 [9] "CCSZ01" "CCSZ01" "BSSZ01" "BSSZ01" "HGSZ03" "HGSZ03" "SESZ04" "SESZ04"
[17] "SLSZ04" "SLSZ04" "PRSZ03" "PRSZ03" "GLSZ05" "GLSZ05" "TMSZ05" "TMSZ05"
[25] "CHSZ02" "CHSZ02"

$ORIGIN_BS
 [1] "11009" "11009" "22501" "22501" "43709" "43709" "47201" "47201" "51071"
[10] "51071" "53041" "53041" "62251" "62251" "67421" "67421" "68091" "68091"
[19] "77329" "77329" "82221" "82221" "96319" "96319" "97079" "97079"

$TRIPS
 [1]  6229  6229  2533  2533  1608  1608 15689 15689 14145 14145 10621 10621
[13]  1874  1874   297   297   982   982   326   326  2871  2871  2996  2996
[25]  1980  1980

$ORIGIN_SZ
 [1] "QTSZ01" "QTSZ01" "JWSZ09" "JWSZ09" "BKSZ07" "BKSZ07" "WDSZ07" "WDSZ07"
 [9] "CCSZ01" "CCSZ01" "BSSZ01" "BSSZ01" "HGSZ03" "HGSZ03" "SESZ04" "SESZ04"
[17] "SLSZ04" "SLSZ04" "PRSZ03" "PRSZ03" "GLSZ05" "GLSZ05" "TMSZ05" "TMSZ05"
[25] "CHSZ02" "CHSZ02"

$ORIGIN_BS
 [1] "11009" "11009" "22501" "22501" "43709" "43709" "47201" "47201" "51071"
[10] "51071" "53041" "53041" "62251" "62251" "67421" "67421" "68091" "68091"
[19] "77329" "77329" "82221" "82221" "96319" "96319" "97079" "97079"

$TRIPS
 [1] 4698 4698 1539 1539  415  415 2204 2204 5163 5163 4001 4001 1431 1431   91
[16]   91  223  223   92   92 2249 2249  164  164  290  290

$ORIGIN_SZ
 [1] "QTSZ01" "QTSZ01" "JWSZ09" "JWSZ09" "BKSZ07" "BKSZ07" "WDSZ07" "WDSZ07"
 [9] "CCSZ01" "CCSZ01" "BSSZ01" "BSSZ01" "HGSZ03" "HGSZ03" "SESZ04" "SESZ04"
[17] "SLSZ04" "SLSZ04" "PRSZ03" "PRSZ03" "GLSZ05" "GLSZ05" "TMSZ05" "TMSZ05"
[25] "CHSZ02" "CHSZ02"

$ORIGIN_BS
 [1] "11009" "11009" "22501" "22501" "43709" "43709" "47201" "47201" "51071"
[10] "51071" "53041" "53041" "62251" "62251" "67421" "67421" "68091" "68091"
[19] "77329" "77329" "82221" "82221" "96319" "96319" "97079" "97079"

$TRIPS
 [1]  6229  6229  2533  2533  1608  1608 15689 15689 14145 14145 10621 10621
[13]  1874  1874   297   297   982   982   326   326  2871  2871  2996  2996
[25]  1980  1980

$ORIGIN_SZ
 [1] "QTSZ01" "QTSZ01" "JWSZ09" "JWSZ09" "BKSZ07" "BKSZ07" "WDSZ07" "WDSZ07"
 [9] "CCSZ01" "CCSZ01" "BSSZ01" "BSSZ01" "HGSZ03" "HGSZ03" "SESZ04" "SESZ04"
[17] "SLSZ04" "SLSZ04" "PRSZ03" "PRSZ03" "GLSZ05" "GLSZ05" "TMSZ05" "TMSZ05"
[25] "CHSZ02" "CHSZ02"
AM_6_9 <- unique(AM_6_9)
PM_5_8 <- unique(PM_5_8)
AM_11_2 <- unique(AM_11_2)
PM_4_7 <- unique(PM_4_7)

Check if duplicates cleared

mpsz_AM_6_9 <- left_join(mpsz,AM_6_9,by = c("SUBZONE_C" = "ORIGIN_SZ"))
mpsz_PM_5_8 <- left_join(mpsz,PM_5_8,by = c("SUBZONE_C" = "ORIGIN_SZ"))
mpsz_AM_11_2 <- left_join(mpsz,AM_11_2,by = c("SUBZONE_C" = "ORIGIN_SZ"))
mpsz_PM_4_7 <- left_join(mpsz,PM_4_7 ,by = c("SUBZONE_C" = "ORIGIN_SZ"))
tmap_mode("plot")
AM_6_9_map <- tm_shape(mpsz_AM_6_9)+
  tm_fill("TRIPS", 
          style = "quantile", 
          palette = "Blues",
          title = "Passenger trips") +
  tm_layout(main.title = "Passenger trips generated Weekday 6am-9am",
            main.title.position = "center",
            main.title.size = 0.7,
            legend.height = 0.45, 
            legend.width = 0.35,
            frame = TRUE) +
  tm_borders(alpha = 0.5) +
  tm_compass(type="8star", size = 1) +
  tm_scale_bar() +
  tm_grid(alpha =0.2) +
  tm_credits("Source: Planning Sub-zone boundary from URA\n and Passenger trips data from LTA", 
             position = c("left", "bottom"))

PM_5_8_map <- tm_shape(mpsz_PM_5_8)+
  tm_fill("TRIPS", 
          style = "quantile", 
          palette = "Greens",
          title = "Passenger trips") +
  tm_layout(main.title = "Passenger trips generated Weekday 5pm-8pm",
            main.title.position = "center",
            main.title.size = 0.7,
            legend.height = 0.45, 
            legend.width = 0.35,
            frame = TRUE) +
  tm_borders(alpha = 0.5) +
  tm_compass(type="8star", size = 1) +
  tm_scale_bar() +
  tm_grid(alpha =0.2) +
  tm_credits("Source: Planning Sub-zone boundary from URA\n and Passenger trips data from LTA", 
             position = c("left", "bottom"))

AM_11_2_map <- tm_shape(mpsz_AM_11_2)+
  tm_fill("TRIPS", 
          style = "quantile", 
          palette = "Reds",
          title = "Passenger trips") +
  tm_layout(main.title = "Passenger trips generated Weekend/holidays 11am-2pm",
            main.title.position = "center",
            main.title.size = 0.7,
            legend.height = 0.45, 
            legend.width = 0.35,
            frame = TRUE) +
  tm_borders(alpha = 0.5) +
  tm_compass(type="8star", size = 1) +
  tm_scale_bar() +
  tm_grid(alpha =0.2) +
  tm_credits("Source: Planning Sub-zone boundary from URA\n and Passenger trips data from LTA", 
             position = c("left", "bottom"))

PM_4_7_map <- tm_shape(mpsz_PM_4_7)+
  tm_fill("TRIPS", 
          style = "quantile", 
          palette = "Greens",
          title = "Passenger trips") +
  tm_layout(main.title = "Passenger trips generated Weekends/holidays 4pm-7pm",
            main.title.position = "center",
            main.title.size = 0.7,
            legend.height = 0.45, 
            legend.width = 0.35,
            frame = TRUE) +
  tm_borders(alpha = 0.5) +
  tm_compass(type="8star", size = 1) +
  tm_scale_bar() +
  tm_grid(alpha =0.2) +
  tm_credits("Source: Planning Sub-zone boundary from URA\n and Passenger trips data from LTA", 
             position = c("left", "bottom"))
#tmap_arrange(AM_6_9_map, PM_5_8_map,AM_11_2_map,PM_4_7_map, asp = 2, ncol = 2)

Transform data to hex friendly

gtest_points = mpsz_AM_6_9 %>%
  # lng/lat value are missing in some records
  filter(!is.na(geometry)) %>%
  st_as_sf(coords = "geometry", crs = 3414, remove = FALSE)

Making grids

hex_grid <- st_make_grid(gtest_points, c(750, 750), what = "polygons", square = FALSE)
grid_sf = st_sf(hex_grid) %>%
  # add grid ID
  mutate(grid_id = 1:length(lengths(hex_grid)))
contents <- st_intersects(hex_grid, gtest_points)

for (i in 1:length(contents)) {
  x = 0
  for (j in unlist(i)){
    x = x + gtest_points$TRIPS[j]
  }
  grid_sf$ridership[i] = x

}
grid_sf$ridership = lengths(st_intersects(hex_grid, gtest_points))
grid_count = filter(grid_sf, ridership > 0)
AM_6_9_map <- tm_shape(grid_count)+
  tm_fill("ridership", 
          style = "quantile", 
          palette = "Blues",
          title = "Passenger trips") +
  tm_layout(main.title = "Passenger trips generated Weekday 6am-9am",
            main.title.position = "center",
            main.title.size = 0.7,
            legend.height = 0.45, 
            legend.width = 0.35,
            frame = TRUE) +
  tm_borders(alpha = 0.5) +
  tm_compass(type="8star", size = 1) +
  tm_scale_bar() +
  tm_grid(alpha =0.2) +
  tm_credits("Source: Planning Sub-zone boundary from URA\n and Passenger trips data from LTA", 
             position = c("left", "bottom"))
AM_6_9_map

tmap_mode("view")

map_honeycomb = tm_shape(grid_count) +
  tm_fill(
    col = "ridership",
    palette = "Blues",
    style = "cont",
    title = "Number of bus stops",
    id = "grid_id",
    showNA = FALSE,
    alpha = 0.6,
    popup.vars = c(
      "Number of bus stops: " = "ridership"
    ),
    popup.format = list(
      ridership = list(format = "f", digits = 0)
    )
  ) +
  tm_borders(col = "grey40", lwd = 0.7)

map_honeycomb
mode(map_honeycomb)
[1] "list"
mode(map_honeycomb)
[1] "list"
mapview_test_points = mapview(gtest_points, cex = 3, alpha = .5, popup = NULL)

mapview_test_points